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BOSE-EINSTEIN CONDENSATION AT FINITE TEMPERATURES: 
MEAN FIELD LAWS WITH PERIODIC MICROSTRUCTURE 


DIONISIOS MARGETIS* 

Abstract. At finite temperatures below the phase transition point, the Bose-Einstein conden¬ 
sation, the macroscopic occupation of a single quantum state by particles of integer spin, is not 
complete. In the language of superfiuid helium, this means that the superfluid coexists with the 
normal fluid. Our goal is to describe this coexistence in trapped, dilute atomic gases with repulsive 
interactions via mean field laws that account for a spatially varying particle interaction strength. 
By starting with the A^-body Hamiltonian, A” ^ 1, we formally derive a system of coupled, nonlin¬ 
ear evolution equations in 3 -h 1 dimensions for the following quantities: (i) the wave function of the 
macroscopically occupied state; and (ii) the single-particle wave functions of thermally excited states. 
For stationary (bound) states and a scattering length with periodic microstructure of subscale e, we 
heuristically extract effective equations of motion via periodic homogenization up to second order in e. 

Key words, quantum dynamics, Bose-Einstein condensation, periodic homogenization, finite 
temperatures, two-scale expansion, mean field limit 
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1. Introduction. The macroscopic occupation of a single quantum state, known 
as Bose-Einstein condensation (BEC)^ was first observed in trapped, dilute atomic 
gases at very low temperatures almost two decades ago mM- This advance has paved 
the way to the precise control of the quantum behavior of ultracold atomic systems; 
for broad reviews on BEC, see, e.g., [^[71|51[TT1[^[iD]. 

From the theoretical viewpoint, three remarkable features of dilute atomic gases 
in current experimental settings are the following, (i) The weak atomic interactions. 
In many experiments, these interactions are short ranged, characterized by an effective 
parameter, the scattering length a ED; for repulsive interactions, a > 0. (ii) The trap¬ 
ping external potential. This is needed in order to keep the atoms together. Because 
of this trap, there is no translation invariance of the Boson system, (iii) Finite but 
“small” temperatures below the phase transition point. The actual temperatures of 
the experiments are very low, in the range of nanodegrees Kelvin Nonetheless, 

the BEC is not complete: particles at the macroscopic quantum state, the condensate, 
coexist with many particles that occupy thermally excited states. In the language of 
superfluid helium, the superfluid and the normal fluid are both present |33j . 

The first feature [(i)], weak particle interactions, enables a systematic theoretical 
treatment of the Boson gas, as this is carried out in the insightful works by Lee, 
Huang and Yang [30l|3T] who modeled BEC in the periodic case; in this setting, 
the condensate is the zero-momentum state and thermally excited states correspond 
to nonzero momenta. In connection to the second feature [(ii)], external potential, 
an ingenious many-body quantum-mechanical treatment of BEC for systems without 
translation symmetry was developed by Wu m- These treatments [5011311113 make 
use of physically transparent approximations for the many-body Hamiltonian. 

In this article, we focus on the the third feature [(in)], finite temperatures, with 
a spatially varying scattering length. Our analysis addresses the following question. 
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What mean field laws are consistent with the many-body quantum dynamics of the 
trapped Boson gas at finite temperatures, sufficiently below the phase transition? In 
particular, we carry out the following tasks. 

• We formally derive, directly from a microscopic model by perturbation of 
nonrelativistic quantum fields, mean field evolution laws for the trapped di¬ 
lute atomic gas at finite temperatures when the scattering length is spatially 
varying. Our analysis combines heuristic techniques adopted from |.?01I811 
1461147) . The emerging description consists of a system of coupled^ nonlinear 
Schrodinger-type partial differential equations (PDEs) in 3-1-1 dimensions for: 
(i) the wave function, $(t, x), of the condensate; and (ii) the one-particle wave 
functions, (l)j{t,x), of the thermally excited states. This system of PDEs has 
spatially-varying coefficients of nonlinear terms. 

• Eor stationary (bound) states, we provide an effective description for this PDE 
system when the atomic interaction has a periodic microstructure of subscale 
e. This e expresses the ratio of the length over which the scattering length, a, 
varies to the particle correlation length. We apply classical homogenization 
theory to semilinear elliptic equations in the spirit of Bensoussan, Lions, and 
Papanicolaou [4], inspired by work of Fibich, Sivan and Weinstein on the 
focusing nonlinear Schrodinger equation |20) . 

In regard to the derivation of the mean field laws, the proof of convergence of 
the high-dimensional microscopic dynamics to low-dimensional PDEs lies beyond our 
present scope. We invoke an uncontrolled yet physically motivated ansatz for the 
many-body Schrodinger state vector. Furthermore, the convergence of the two-scale 
expansion in the periodic homogenization program is not pursued here. Approximate 
solutions of the PDE system are intended to be the subject of future work. 

1.1. Motivation. Our work is motivated by experimental advances in control¬ 
ling properties of ultracold, trapped atomic systems. Evidently, Bose-Einstein con¬ 
densates form systems of considerable promise for precision metrology which may offer 
improvement over conventional interferometry [131125112^ . Atomic gases undergoing 
BEG can also be used to probe or emulate properties of condensed-matter systems; 
e.g., the presence of quantum phase transitions such as the transition between a su¬ 
perfluid and a Mott insulator [36]. In these applications, the short-ranged atomic 
interactions are known to play an prominent role. These interactions can be tuned 
externally, e.g., by optical means near a Feshbach resonance [9l[T2l|45]. 

These ongoing experimental advances suggest the need to study effects intimately 
connected to the atomic interactions in non-translation-invariant settings. A key the¬ 
oretical aspect of this problem is the possible effect on the description of the Bose- 
Einstein condensate of a spatially varying scattering length, a. This possibility is the 
subject of the formal analysis in [34] . in which mean-field and beyond-mean-field im¬ 
plications for the condensate are studied at zero temperature. The present article is 
meant to be a first step toward extending the perturbative analysis of [34] to the case 
with finite temperatures, when thermally excited states influence the condensate. 

1.2. Overview: Mean field limit at zero temperature. It is useful to recall 
that a system of N interacting Bosons evolving at zero temperature is described by a 
symmetric wave function obeying the iV-body Schrodinger equation. This description 
requires the use of 3iV spatial coordinates and, thus, is deemed as unwieldy for physical 
prediction if iV ^ 1. It is desirable to reduce the many-body evolution to PDEs for 
variables defined in lower dimensions. A related variable is the condensate wave 
function, which typically lives in 3 spatial dimensions. 
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For repulsive atomic interactions, this reduction usually leads to a mean field 
law: a defocusing, cubic nonlinear Schrodinger equation, herein referred to as the 
Gross-Pitaevskii-Wu equation (GPWE) [231[3nillZ]- In this PDF, the coefficient of 
the nonlinear term is proportional to the scattering length, a; cf. section [TM It has 
been rigorously shown by Elgart, Erdos, Schlein and Yau [TtHII] how the reduced, 
one-particle evolution, including the scattering length, properly emerges as fV —>■ oo 
from the many-body dynamics by use of the Bogoliubov-Born-Green-Kirkwood-Yvon 
(BBGKY) hierarchy for particle marginal densities via Spohn’s formalism [U] . 

This mean held limit holds at extremely low temperatures, if particles of the nor¬ 
mal component (“thermal cloud”) of the atomic gas on average have a negligible effect 
on evolution. As the temperature increases, while remaining well below the phase 
transition point, many particles can occupy thermally excited states. The GPWE for 
the condensate must be modihed to account for the occupation of such states. 

1.3. On past works. The stationary case of hnite temperatures well below 
the phase transition was treated systematically by Lee and Yang in the periodic set¬ 
ting [CT32] . Specihcally, the fraction, of particles at the condensate was introduced 
as a hxed parameter in an approximation scheme for the A^-body Hamiltonian; this ^ 
was hnally determined via a statistical average as a function of temperature, T [S]. A 
key idea in the approximation scheme was to treat the fluctuations about the average 
number, A^(l — ^), of particles out of the condensate as small. This idea permeates 
our present treatment; see also [Mj . 

In the early 1980s, studies of Bosons in non-translation-invariant settings were 
motivated by experiments involving spin-polarized atomic hydrogen inside magnetic 
traps. These studies focused on bound states. For example, by minimizing a certain 
temperature-dependent free-energy functional of the self-consistent Hartree-Fock and 
Bogoliubov theory, as described, e.g., in m, Goldman, Silvera, and Leggett [52] 
obtained equations of motion for the time-independent condensate and single-particle 
excitation wave functions. A similar method of variational character was followed by 
Huse and Siggia [28|, who systematically derived equations of motion that retain the 
orthogonality between the condensate and thermally excited states. 

Another approach in stationary settings by Oliva m invoked the local-density 
approximation, by which the dependence of requisite functionals on the local particle 
density was assumed to be that of the homogeneous system; see also [1^. Later 
works, e.g.. [21] . made use of Popov’s diagrammatic approach [41] for the homogeneous 
gas; cf. [40] for a review. More recent treatments rely on PDE models that couple 
the condensate with the heat bath via non-equilibrium methods of quantum kinetic 
theory [51[^ [T51H5] . For a review of various models, see, e.g. [22]. 

1.4. Main approach and mean field evolution laws. Our approach to de¬ 
riving the mean field limit relies on: (a) a perturbative approximation scheme for the 
A^-body Hamiltonian in the spirit of [31] in the non-translation-invariant framework 
of [47] for zero temperature; and (b) a time-dependent ansatz for the A-body state 
vector that incorporates the distribution of particles over thermally excited states 
according to Bose statistics. The scheme is expected to hold when: 6 = \/qo? <C 1, 
where q is the gas local density and a is the scattering length; and fluctuations about 
the average number of particles at each thermally excited state are small enough m 
(see section [331 for details). The (small) parameter 5 expresses the gas diluteness. 

Our program is carried out in the Bosonic Fock space F (defined as the Hilbert 
space with an indefinite number of Bosons), by making use of the fraction, of the 
particles at the condensate as a fixed parameter in the approximation scheme [31] : 
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cf. |35j . Our treatment yields low-dimensional evolution laws; and, thus, is in principle 
distinct from previous treatments of bound states such as the variational approach 
of [^1^ , the local-density approximation of [TOIISZ] : and the finite-temperature the¬ 
ory of [21] . 

By starting with the A^-body dynamics, we derive the following evolution PDEs for 
the condensate wave function, $(t,a;), and single-particle excitation wave functions, 
(j = 1, 2, ...): 

(1.1a) idt^ = {-Aj, -h Ve(x) + g{x)[Qs{t,x) + 2Qnlt,x)] - , 

idt(i)j = {-A,;, -t Ve{x) + 2g{x)[Qs{t,x) + Qn{t,x)] - 
(1.1b) - 6j(t)$(t,x) ; 

which for bound states reduce to the corresponding laws in [^. We assume that 
ll®(^G)lli2(R3) = N, ((^i(t, •),4 (t, •))i2(R3) = 5ij and ($(t, •), •))l 2(R3) = 0. In 

the above, A^, denotes the 1-particle Laplacian; Ve{x) is the external potential, as¬ 
sumed to be smooth and grow rapidly for large |a;|; g(x) = 87ra(x) where a(x) is 
the spatially varying scattering length; is the mean field condensate (superfluid) 
density, Cl'hp; Pn is the thermal cloud (normal fluid) density [cf. (14.101) ]: bj(t) comes 
from the orthogonality of each (j)j to $ [cf. (14.361) ] ; and 

(1.2) dxg{x)\^{t,x)\'^ . 

Jr3 

The parameter ^ can in principle be determined as a function of temperature, T, via 
imposing Bose statistics for the occupation numbers of particles at thermally excited 
states (see section[5]) [22]. At zero temperature, when ^ = 1 and pn = 0, (ll.la|) readily 
reduces to the GPWE [^15^1 ITT] . 

For bound states, we set ^{t,x) = e“*^*$(a:) and ^j{t,x) = e~'‘^^*(j)j{x) where E 
and Ej are the associated energies per particle (section [5]) . 

1.5. Periodic microstructure and homogenization. Following [^, we in¬ 
troduce the function 

(1.3) g{x) = 87ra(x) = go[l + A{x/e)\ , go = Sttuo > 0 , 

in (II. 1[) : A{y) is a smooth periodic function of zero mean and period equal to unity. If 
I sc is the length over which the scattering length varies and Ic is the particle correlation 
length, then e = Isc/lc- We assume that: Igc is small compared to other length scales 
of the problem such as Ic, a, the thermal de Broglie wavelength and mean interparticle 
distance; and the external potential, Ve{x), does not depend on the subscale e. 

For bound states, we formally apply the two-scale expansions 

OO 

(1.4) u{x) = u^°\x, x) +u^'^\x, x) , x = xle, 

k=l 

where u = ^ or pj] = 0(1) as e 0. Then, we reduce (II.1|) to a set of e- 
independent equations of motion for u^^\x,x) with fc = 0, 1, 2 (section [6]) [4ll^. 

This microstructure has been inspired by a model of Fibich, Sivan and Wein¬ 
stein |2^. Notably, our formal analysis pays attention to system dm of defocusing 
nonlinear Schrbdinger-type PDEs, to be contrasted with the single focusing PDF 
studied in [2D] . 
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1.6. Limitations. The present work points to several pending issues. Our ap¬ 
proach is heuristic, placing emphasis on modeling aspects of the Boson gas; the rigor¬ 
ous analyses of the mean field limit and homogenization procedure are not pursued. 
We do not carry out any numerical simulation for the derived PDEs; this task is 
left for future work. In the periodic homogenization program, we focus on stationary 
states; the time-dependent problem, in which two time scales should be accounted for, 
remains unresolved. The rapid spatial variation of the scattering length, a, can be 
driven by a corresponding variation of the external potential, 14 , in optical traps [3]. 
The treatment of the microstructure of both 14 and a is the subject of near-future 
work. We herein restrict attention to mean field laws, leaving out effects of quan¬ 
tum fluctuations due to the excitation of atoms in pairs from the condensate m- 
The description of this process for a scattering length of periodic microstructure at 
zero temperature is provided in [34) . In the present case of finite temperatures, it is 
expected that PDEs (El) would need to be coupled with a nonlocal PDE for the req¬ 
uisite “pair-excitation kernel” which should depend on the parameter 4 Tbe modeling 
of this more demanding process is not touched upon in our treatment. 

1.7. Structure of article. Section 2 presents some terminology and notation. 
In section [3l we review concepts of the second quantization, which are needed in our 
derivations. Section 0] focuses on the derivation of mean field evolution laws. In sec¬ 
tion [SJ we focus on the formalism of stationary states. In section (B) we homogenize 
the stationary equations of motion for a scattering length of periodic microstructure. 
Section [7] concludes our article with a summary of results and open problems. 

2. Notation and terminology. We adopt the following conventions through¬ 
out this article, in the spirit of [34) . 

• C is the complex plane, Z is the set of all integers, N = {0, 1, ...}, and 
N* = N \ {0}. The star (*) as a superscript implies Hermitian conjugation. 

• is the space of symmetric (Bosonic) functions on which are 
invariant under permutations of the particle spatial coordinates, (cci,..., Xn)- 

• (F, G) is the one-particle L^-inner product, F(x)G{x) da:, with induced 
norm ||F||. The scalar product in the Fock space, F, is denoted (•,-)f (see 
section [321); the induced norm is || • ||f. 

• One-particle integrals with unspecified integration ranges are meant to be 
integrals on Integration by parts yields vanishing boundary terms at oo. 

• denotes the d-dimensional unit torus (cell), where d = 3 for our purposes. 
Functions that satisfy A{x -\- e^) = A{x) for all x = (x^, x'^, x^) G R^ and 
A: = 1, 2, 3, where e^’s are unit Cartesian vectors, are referred to as 1-periodic. 
(4.) is the cell average (or mean) of the 1-periodic A{x). 

• denotes the Sobolev space for fc = 1 and p = 2, with dual space 

and id^v is tbe space of 1 -periodic functions with zero average. 

• The dual space = {/ g id“^(T'^)| (/) = 0} is the Hilbert space 

equipped with the inner product (/, h)= 

denotes the ida;j,^-norm of the 1-periodic A{x); typically, d = 3. 

• Let (the 1-periodic) F be in L^(T^). If F has zero mean, then {—Ax)~^F{x) = 

s > 0, where ^ i® the Fourier series 

representing F. Note that ((—A)“'*F) = 0, d = (d^, £^) G Z®. Here, £ ■ x 

is the (Euclidean) scalar product of the 3-dimensional vectors £ and x. 

• Writing f = 0{g) [f = 0 ( 5 )], where / and g are functions, means that \ f/g\ 
is bounded by a nonzero constant (tends to zero) in a prescribed limit. / ~ 5 
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is used loosely to imply f — g = o{g). If A and B are operators on F, ^ ^ S 
means that ('I',^'h)][r ^ ('i',i3'h)F where ^(t) is any Schrddinger state vector 
with components on n-particle sectors of F and n ^ 1. 

• is the k-th order coefficient in the two-scale expansion for u = (j)j; 

cf. (II.4|) . We use the symbols and fj to denote the slowly-varying parts 

of and (section analogous notation is used for the coefficients, 

and in the expansions for the superfluid and normal fluid densities. 

• The dot on top of a symbol denotes time derivative. 

3. Many-body Boson dynamics: Background. In this section, we describe 
the many-particle Hamiltonian (section lXTl) . review basics of the Fock space formalism 
(sectioning, and outline the idea of the many-body perturbation scheme (section ing : 
for related reviews, see [3113]. 

3.1. Many-particle Hamiltonian. The starting point is the Hamiltonian op¬ 
erator, H]s[, of N Bosons [Hn : Lg(]R^^)], viz., 

N 1 

(3.1) Hi\j = 'y ( [—Aj -|- Ve{xj)] — 'y ( Vint {xj , Xj) {Xj G M ) , 

1=1 

which we will refer to as the “PDF Hamiltonian”. The units are chosen such that 
h = 2m = 1 {h: Planck’s constant, to: atomic mass) and xj are particle positions in 
The interaction potential Vint is positive, symmetric, and compactly supported. 
We assume that the external potential, Ve(x), is positive, smooth and increasing 
rapidly for large |a;| with Ve{x) —>■ oo as |a:| —>■ oo. 

In the following, we replace Vint by 

(3.2) Vintixi^Xj) ^ V{xi,Xj) = g{x^) S{xi - Xj) , 

where g{x) = 87ra(x) > 0 and 6{x) is the Dirac mass in Hence, the scattering 
length, a{x), enters our model as an ad hoc parameter in the spirit of [47) . This is to 
be contrasted with the rigorous approach in [IZHIS] where the (constant) scattering 
length emerges from the limiting procedure. 

3.2. Fock space formalism. The Bosonic Fock space, F, consists of vectors v 

formed by sequences of n-particle symmetric wave functions; G Lg(R^"), 

n G N [ 3 . In particular, D := {1, 0, ...}, D G F, is the vacuum state, which has 
no particles. The A^-particle Schrddinger state vector on F is represented by ^(t) = 

with = ^N5n,N', ^Ar(-,t) G Ls(R^^). The scalar product on F is defined 
by (ui, V 2 )v = 

For the one-particle wave function / G L^(R^), the creation and annihilation 
operators and a/ on F are defined by 

n 

{a*fvY'^\xn) = /(Xj)t;^”“^^(a:i,.. .,Xj-i,Xj+i, ... ,x„) , 

i=i 

(a/u)(”)(x„) = y/n-\-l y dx/*(x)u*^”+^)(x,x„) ; 

Xn = (xi,... ,x„). Thus, [«/,«*] := a/a* - a*gaf = {f,g) and [af,ag] = 0. The 
operators and a/ create and annihilate a particle at state /, respectively. The 
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Boson field operator 4’{x) and its adjoint, are defined via 

(3.3) ^ V’(a;) ; 

ip*(x) creates (annihilates) a particle at position x and is time-independent in 

the Schrddinger picture. Note the canonical commutation relations ['lp{x),^p*{'y)] = 
5{x—y) and [V’(x),')/’(?/)] = 0. Evidently, 'ip{x)n = 0. The (Hermitian) particle number 
operator. A/”, on F is defined by A/” = / dxilj*{x) ipix) which obeys (A/'u)^") = nv^^\ 
for every u G F and n G N. 

In view of this formalism, the PDE Hamiltonian, Hn, corresponds to the operator 
7^ : F —^ F where specifically, by (13.11) and (13.21) . 

(3.4) 7^ = J dx'ip*{x)[-A^+ Ve{x)]i^{x) J dx'ip*{x)^g{x)i){x)^ . 

We will henceforth use this TL in the place of the PDE Hamiltonian. 

3.3. Many-body perturbation. The starting point of the perturbation scheme 
is to split tp{x) as [47] 

(3.5) -tpix) = -ipoit, x) + x) , 

where is the Boson field annihilation operator for the condensate and tp±{t,x) 
is the Boson field annihilation operator in the space orthogonal to the condensate; 
J dx^*{t,x)ip±{t,x) = 0 where is the condensate wave function. By defining 
aoit) = with ||$||| 2 (r 3 ) = N, we set 

(3.6) ^|;o{t,x) = N~^^'^ao{t)^it,x) ; 

[ao(t), 09 ( 1 )] = 1 and ao(t)D = 0. The operator 'ip± can be expanded as 

00 

(3.7) tp±{t,x) ='^(l)j{t,x)aj{t) , 

i=i 

where (slightly abusing notation) aj{t) is used in place of a^.{t) to denote the an¬ 
nihilation operator for a particle at state (j)j- Together with its adjoint, every 

aj{t) satisfies the canonical commutation relations for Bosons: [a*{t),aj{t)] = Sij and 
[ai{t), aj{t)] = 0. Hence, the operator a*{t)aj{t) on Fock space represents the number 
of particles at state . 

The perturbation analysis relies on expanding % in powers of and V'l. This 
expansion needs to be combined with the following formal statement. 

(3.8) ao(t)ao(t)-f J dx ^p]_{t,x)'tpj_{t,x) = N , 

which sets the particle number operator, M = J dx'ij)*{x) ipix), equal to N (times the 
unity operator). This statement poses the constraint that the many-body Schrddinger 
state vector, 4'(t), be an eigenvector of Af with eigenvalue equal to N for every f > 0; 
A/’4'(t) = In (|3.8p . J dx'ijj]_{t,x)'ip±{t,x) is the operator for the number of 

particles out of the condensate. 

At this stage, it is useful to provide the following definition. 
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Definition 3.1. (Condensate fraction, The condensate fraction, f,, is defined 
as the parameter equal to {^{t), {aQ(t)ao{t)/{t))f, where 'l'(i) is the Schrodinger 
state vector, G F; for r>0, 0<^<1 (see section 0- 

This ^ will be invoked as a T-dependent parameter in the many-body perturbation 
scheme of section 01 By (13.81) . the average total number of particles at thermally 
excited states is ('k(t),/ dxf:]_{t,x)tp±{t,x)^{t))w = N{l — ^). 

Definition 3.2. (Quantum fluctuations in thermally excited states). Consider 
the operator-valued occupation number Mj{t) = a*(t)aj(t) (j G N*j on F for each 
thermally excited state (fj. The operator 5Mj(t) := {'^{t),Mj'^{t))r—Mj onF describes 
quantum fluctuations in the number of particles at state 4>j. In this vein, the operator 

5Mx{t) '■= .^(1 ~ ~ J x) 

expresses quantum fluctuations in the total number of particles occupying thermally 
excited states. 

At this point, we introduce the following assumption, which serves as a rule for 
our heuristics. 

Remark 3.3. The number of particles at each thermally excited state (j>j undergoes 
small fluctuations about its average value. In other words, the 6Afj(t) of Definition \3.‘A 
is treated as “small”, in some appropriate sense, in our approximation scheme (see 
section^for details). 

A consequence of the above assumption is the following statement. 

Remark 3.4. The operator 6Af± of Definition \3.‘A is treated as sufficiently “.small” 
in our approximation scheme; by (13.81) . ditto for the operator — aQ{t)ao{t). 

At T = 0, where ^ — 1, the mean field limit is consistent with treating as 

small, in some appropriate sense, and retaining up to linear-in-^_L, ip*) terms in R via 
a tensor-product ansatz for di |341l47j . This program must be modified for nonzero 
temperature, T f 0: Because of the presumed presence of many, up to 0{N), particles 
at the thermal cloud, quadratic-in-ip±-and-f')_ terms must appropriately be retained 
in TL (see section HU) . These terms will be collected in view of Remarks 13.31 and 13.41 

4. Derivation of mean field evolution laws. In this section, we derive mean 
field laws for the condensate and one-particle thermally excited states for given f and 
occupation numbers, Uj, of such excited states. For this purpose, we replace the many- 
body Hamiltonian "H on F by an operator that contains up to certain quadratic-in-'!/)_L- 
and-ip*^ terms, consistent with the presence of iV(l — particles at thermally excited 
states, where 0 < ^ < 1 and N 1 (section0T]) . This replacement involves a sequence 
of steps enforcing constraint (13.81) under R.emarks l3.3l and l3.41 Evolution PDEs (11.11) in 
3-1-1 dimensions are extracted by enforcement of the many-body Schrodinger equation 
on an uncontrolled but self-consistent ansatz for the Schrodinger state vector 'l'(t) 
(section 0^. 

4.1. Approximations for iV-body Hamiltonian. The Fock space Hamilto¬ 
nian (13.4|) is recast to the form [35] 

n=T+v, 


(4.1) 
where 

(4.2) 


T= / dx{\V'fi{x)(^ -^Ve{x)ip*{x)il}{x)} 
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and 

(4.3) ^ ^ ^ J dx Ip* {xfg{x)ip{xf . 

4.1.1. Interaction operator, V. First, we concentrate on the interaction op¬ 
erator given by ()4.3|) . By use of splitting ()3.5p , V is directly decomposed as 

(4.4) V = Vo -b Vi -t V2 -t V3 + V4 , 


where 'tp]_{t,x) and ip±{t,x) explicitly appear k times in Vk {k = 0, 1, 2, 3, 4). Specif¬ 
ically, we obtain the terms 

(4.5a) Vo = , 

Vi =N~^^‘^ao(t) J dxg{x)\^{t,x)\‘^{^*{t,x)aQ{t)'ipj_{t,x) 

(4.5b) + ^{t,x)tp*j_{t,x)aoit)}ao{t) , 

V 2 = J dxg{x){^*{t,x)‘^aQ{t)‘^-ipj_{t,x)^ + ^{t,x)‘^pj*j_(t,xfao(t)‘^ 

(4.5c) -b 4|l>(t, x)f aQ{t)ip*j_{t, x)ip±{t, x)ao{t)}, 

V 3 = J dxg{x){^*{t,x)al{t)ip*j_{t,x)ip±{t,x)‘^ 

(4.5d) + ^{t,x)'ip*j_{t,xytp±{t,x)aoit)} , 

(4.5e) V 4 = 5 y dxg{x)^p*J_(t,xftp±(t,xf , 

where ({t) is defined by (II.2L Equations (14.51) need to be further manipulated for 
the extraction of their full dependence on ip± and ip*^. Specifically, constraint (I3.8p 
imposes a relation between OgOo and that has not yet been accounted for 

in (113]) : this task will be carried out in the sequence of steps described below. 

Next, we show that all terms Vk contribute to the mean field limit, in contrast to 
the zero-temperature case [34]. This has been ascertained in the periodic setting at 
T > 0 in which thermally excited states correspond to nonzero momenta m- 

The Vi of (|4.5bp is the simplest term to deal with since in this term we can simply 
set ag(t)ao(t) = iV^, according to Definition 13.11 and Remark 1X11 Thus, we make the 
replacement 

Vi - J dxg{x)\^{t,x)f{^*{t,x)al{t)ip±{t,x) 

(4.6) +^{t,x)ip*j_(t,x)aoit)} . 


For Vo, the procedure of approximation takes into account (13.81) and Remark 13.4 

Vo = iA^“^C(0{(ao(Oao(i))^ - aS(Oao(i)} 

^ ^N-^Cit){a*oit)ao{t)r 




Nil-0- J dx'ip*j_it,x)'ip±it,x) 




Nil-0- J dxip’^it,x)'p±it,x) 


= hN-\it) 


N^^i2-0-2N^ J dxr±it,x)i;±it,x) 


(4.7) 
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The term aQ{t)ao{t) inside the braces in the first line of (14.71) contributes 0{N^) and 
0(1) terms which can both be neglected in comparison to the bracketed terms of the 
last line. (We treat 1V(1 — — J dx as small compared to N^, by Remark 15^ 1 

Next, we address the term V 4 of (I4.5e|) . In view of expansion (13.71) . this term is 
rewritten as 

00 

V 4 = i a*(t)a*(t)ak(t)ai(t) 


X j dxg(x)(j)* (t, x)^* (t, x)(j)k(t, x)^i (t, x) 



recall that Mj{t) = a*(t)aj(t) is the operator for the number of particles at state (j}j. 
The approximations in (14.81) rely on: (i) the hypothesis that the major contribution 
to the sum of the first line comes from terms with (fc, 1) = {i,j) or (j, i); and (ii) the 
identity a*a*aiaj = (a*ai)(a*aj) — Sija*ai = JViJVj — 5ija*ai and neglect of 5ija*ai 
because of its relatively small contribution. Hypothesis (i) is responsible for the 
appearance of the extra factor of 2 in the third line of (14.8p . This hypothesis leads 
to the normal fluid density, as an emergent effective potential in the equations of 
motion [see (14.101) ]. in the spirit of the mean field approach that we have adopted. 

To further simplify V 4 , we apply Remark [331 As discussed in [31], it is desirable 
to start with an Wbody state vector, 'k(t), with the occupation numbers 

(4.9a) no = {^(t), {a5(t)a5(t)}^(t)) = , n, = {^(t),Nj(t)^(t)) = 0(1) , 

for the states $ and respectively, where 

CXD 

(4.9b) ^ n, = iV(l - 0 ; 

and treat Afj — nj as small compared to nj in the approximation scheme. In section l431 
we identify nj with occupation numbers at equilibrium. Thus, (14.81) becomes 






i=i 


i=i 


dxg(x) 

dxg(x) 


i=i 

00 

i=i 


-Y +^Y l4(^>2;)pAf)-(t) 


1=1 


1=1 


-Y +2tlj*j_(t,x)tlj±{t,x) 

1=1 
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By (13.71) . in the last line we invoked the relation 

OO 

OO 

At this stage, it is natural to define the mean field density for the normal fluid as 

OO 

(4.10) gn{t,x) \<j)j{t,x)\'^nj . 

1=1 


Thus, Vi is conveniently written as 

(4.11) Vi'^ J dxg{x) gn{t,x)[-gn{t,x)+ 2ij;*j_{t,x)ip±{t,x)] . 

Next, we address the term V 2 of (I4.5cl) . This term is split as 

(4.12) V2=V21+V22, 

where V 21 describes pair excitation, i.e., the scattering of particles in pairs from the 
condensate to other states and vice versa, viz., 

V 21 = J dxg{x){^*{t,xfaQ{t)‘^'>p±(t,xf 

(4.13) +^{t,xf'ilj*j_{t,xfaoitf}, 
while 

(4.14) V 22 = 21V“^y dxg{x)\^(t,x)\‘^aQ{t)'>p*j_{t,x)'>p±{t,x)ao{t) . 

In the following treatment, we neglect V 21 , since its retainment has been shown to 
lead to beyond-mean-field effects at T — 0 [23l|47]. In contrast, V 22 contributes to 
the mean field limit, and must be retained and simplified by (13.81) under Remarks 13.41 
and 13.31 


V 22 = 2A1-1 + 


J dxi’*j_{t,x)ipj_{t,x) 


X / dxg{x)\^(t,x)\‘^{gn{t,x) + [ip*j_{t,x)'ip±{t,x) - p„(t,a;)]} 


(4.15) 


2N- |7V(1-^) J da; 5 f(a;)|d>(t, j;)| p„(t,a;) 

+ A^Cy <^xg{x)\^{t,x)f'il;*j_{t,x)'iljj_{t,x) 

^y9iy)\^(t,y)\'^gn{t,y) / dx'ii!*j_{t,x)ip±{t,x)}. 
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A similar procedure can be applied to the V 3 of (I4.5dl) : 


V 3 = dis(i) 


(j)*{t, x)(i)j {t, x)^k (t, x)a* {t)aj {t)ak (t) 


i,3,k 


x) X! {t, x)a* {t)a*{t)ak (t)ao {t) 


i,3,k 


J dxgix) 


{t, x)aQ (t) ^ (j)* {t, x)^, {t, x)^j {t, x)a* {t)ai{t)aj (t) 




+ x) ^ (t, x)(f)* (t, x)(j)j (t, x)a* {t)a* (t)aj (t)ao (t) 




= 2N J dx g{x) 


$*(t,a:)ao(t) ^ \(j)i{t,x)\‘^Af,{t)ipj_{t,x) 


i=l 


^{t,x)ilj*j_{t,x)ao{t) ^ |(^i(t,x)pA/'i(t) 


i=l 


- 2iV dxgix) gn{t,x)[^*{t,x)al{t)'ilJx{t,x) 

(4.16) + l>(t,x)V'l(t,x)ao(t)] . 


Therefore, the desired approximation to the V of (14.31) is given by (14.41) and (I4.12p 
together with expressions (14.71) . (14.61) . (14.151) . (I4.16p . (14.111) . and V 21 — 0: 

(1“0^]+2(1-C) J dxg{x)\^it,x)\'^gn{t,x) 

- J dxgix) £i„(t,x)^| 

+ J dx g(x){C|$(t, x)p + 2gnit, x)}{<&*(t, x)ao(t)'0x(t, x) 

+ $(t,x)V’l(t,x)ao(t)} 

+ J da: |-CC(i) + 25(x)[^|$(t,x)p + p„(t,x)] 

(4.17) -2N-^J dygiy)\^it,y)\'^gnit,y)'^ijj*j_it,x)ij;j_it,x). 

4.1.2. Kinetic energy operator, T, and resulting Hamiltonian. The ki¬ 
netic energy operator, T, of ()4.2p is explicitly written as 

T = N[CAit) + Ceit)] + N~^/^aoit) J dx [-A^it,x) + Veix)^it,x)]-ip*j_it,x) 

+ JV-'/=.;(,) j ,) + VMi -«. «. X) 

(4.18) -f J dx{\V'lp±it,x)\^ + [Veix) - (Ait) - Ceit)]lp*j_it,x)'lp±it,x)} , 
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where C{t) is defined by (11.21) and 

(4.19) CAit)=N-^ J dx|V$(t,x)|2 , Ce{t)=N-^ J dxVeix)mt,x)\^ . 

Accordingly, the iV-body Hamiltonian is approximated as 

(4.20) -H --H/=-Ho + ■Hi +-^2 , 
where 

Ho = N[U{t) + Ce{t)] + ^Ncm - (1 - o'] 

(4.21) +2(1-0 J dxg{x)\^{t,x)\'^gn{t,x) - J dxg{x) gn{t,x)'^ 

is a c-number (multiplied by the unity operator), while the operators Hi and H 2 
depend on '0j_, 

'Hi = J da; {—A$(t, x) + 14(a;)$(t, x) 

+ 5(x)[0$(t,x)p + 2 gn{t,x)]^{t,x)}ip*j_{t,x) 

_l_ 7V“^/^aQ(t) J dx {—A<l>*(t, x) + t4(x)$*(t, x) 

(4.22) + g{x)[^\^{t,x)f + 2gnit,x)]^*{t,x)}'4’±(t,x) , 


(4.23) 


H 2 = y dx ||Vi/;_L(i,a;)P + 

+ 2g(x)[C|$(t,x)p + £>„(t,x)] 

- 2N-^ [ dyg{y)\^{t,y)\'^gn{t,y) 


Ve{x) - Uit) - Ce{t) - ^at) 

'lp*j_{t,x)'lp±{t,x) 


We view the resulting H/ as the “lowest-order perturbation” for Hamiltonian H in 
the present context of finite temperatures. 

4.2. Mean field PDEs. Next, we combine (I4.20ll - (l4.23p with an ansatz for 
the Schrodinger state vector ^(t) in order to derive PDEs (11.11) . For zero temper¬ 
ature, T = 0, the mean-field approximation for d’(t) is the tensor product 'l'(t) ~ 
(iV!)“^/^a5(t)^D, in which all particles occupy the one-particle state $. For the 
present case with finite temperatures, we modify this ansatz to 


(4.24) 


^(t) 






n 




K-!)V2 


n 


where 9{t) is a global phase to be determined and {nj}^i amounts to the occupa¬ 
tion numbers of the thermally excited states; ao(t) and aj{t) are described in (13.61) 
and (I3A1) . We consider ansatz (14.241) as a mean field approximation consistent with 
the lowest-order H/ of (j4.20|) - (j4.23|) . for small fluctuations of the j-th state occupation 
numbers. A/) , about their averages. The enforcement of this consistency through the 
many-body dynamics of the Schrodinger equation should yield the desired evolution 
laws for $ and 
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Several comments on (14.241) are in order. 

(i) Recall that rij satisfy (I4.9bl) and nj G N. Therefore, only a finite number of 
these Tij can be nonzero. Thus, product (14.241) is finite. 

(ii) Since the nj are integers, dnj/dt can only be delta functions. However, no 
such delta functions are present in the Hamiltonian. Hence, in the present scheme 
all rij must be time-independent. We henceforth set rij = to emphasize the 
equilibrium-like character of these occupation numbers. 

(iii) This time independence of is related to a point discussed in m for zero 
temperature, namely, that the underlying approach is valid only at a moderate time 
scale. This restriction holds for nonzero temperatures as well. The issue of shorter 
and longer time scales remains open. 

(iv) The factor is introduced for later algebraic convenience. By suitable 

choice of the phase, 9[t), the equations of motion for $(t, x) and (j)j{t,x) can be 
relatively simplified. 

To this lowest order of perturbation theory, by neglect of pair excitation, the 
many-particle Schrodinger equation reads 

(4.25) , 


where TLi is approximation (I4.20|) for the Fock space Hamiltonian, TL. To facilitate 
the operator algebra implied by (14.251) . define the following state vectors. 


(4.26a) 


^(o)(i) 




[(Ar^)!]i/2 n 


eq 

(nf!)i/2 


n 


and, for f G N*, 


(4.26b) 






0 ^ TT ' Q 

mv-v/^ * (n0!)i/2 l^i(nf!)V2 ’ 


Clearly, 'I'(i)(t) = 0 if = 0. Thus, there is only a finite number of nonzero 'I'(i)(t) 
by this formalism. The left-hand side of (14.251) equals 

CXD 

= d(t)^(t) -hfdo(t)4'(o)(t) -b^ fd*(t)4'(j)(t) . 

3 = 1 


By virtue of (I4.21I) - (I4.23I) and after some algebra, (14.251) reads 


(4.27) emit) + + Y. = > 


where 


B^it) = ia^it) — N j dx { — A$(t, x)-b t4(x)$(t, x) 

+ C5(a;)l^(C x) + 2g{x)gnit, x)$(t, x)} V’KC x) 


(4.28) 
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and, for j G N*, 

J dx { - A^*{t,x)+ Veix)^*{t,x) 

+ ^g{x)\^{t, x) {t, x) + 2g{x)gnit, {t, x) (f, x) 
- y dx |-A4 (f,x) + Ve{x) - Ca(<) - Ce{t) 
+2^g{x)\^{t,x)f +2g{x)Qn{t,x) 

(4.29) - 2N-^ [ dyg{y)\^{t,y)\‘^gn{t,y) (l)j{t,y)\'ilj*j_{t,x) . 


By (I4.27I) - (I4.29L the operators Bj{t) should satisfy 
(4.30) B*o it) = CO (<)aS (t) , B* (t) = c, {t)a* (t) , 

where co(t) and Cj{t) are such that 


(4.31) 


0{t) + N^co{t) + ^ n^Cj{t) = Ho ■ 
1=1 


Equation (14.311) is the only condition that 9{t), co{t), and Cj(t) must satisfy. 

Equations (14.301) with (14.311) yield PDEs for <l>(t, x) and ^/{t, x), in view of (I4.28P 
and (14.291) . Indeed, taking into account that ao(t) = J dx^{t,x)ip*{x) and 

^l}]_{t,x) = V'*(x) — lV“^/^ao(t)‘h*(t, x), we write relation (14.301) for yBo(t) as 

+ J ‘^y IHx - y) - N~^^{t,x)^*{t,y)] 

x[- Ay + Veiy) + ^giy)\^{t, y)P + 2g{y)gn(t, y)] l>(t, y)| = 0 . 

Thus, the quantity in the braces of the last expression must vanish, leading to 
idt^{t,x) = {-Aj, + V;(x) + 5 (x)^|$(<,x)p + 2g{x)gn{t,x) 

— —Co(t) + CA(t) + Ce(i) + CC(^) 

(4.32a) + 2N-^ J dyg(y) £>„(t, y)|$(<, j/)p |$(t,x) . 

In the same vein, by (14.291) and (I4.30p for B*{t) (j G N*), we obtain the PDE 


idt^j(t,x) = 


-A^ + Ve{x) + 2g{x)^\^{t, x)p + 2g{x)gnit, x) 

- CA(t) - Ce{t) - CC(i) + Cjit) 

-2A"i / dyg{y) gn{t,y)\d>{t,y)\'^ 4(^,x) 


(4.32b) 


-N ^d>{t,x) / dy^*{t,y)g{y)^mt,y)\^^j{t,y) , 
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where the last term in (I4.32bl) results from the orthogonality of $ and . 

Equations (14.321) can be further simplified. By (|4.3ip . we write 

OO 

9{t) + N^co{t) + ^ 

1=1 

= iV[CA w + CeW] + - (1 - en 

+ 2 ( 1-0 j dxg{x)\<^{t,x)\'^Qn{t,x) - J dxg{x)Quit, xf . 

By comparing the last equation to PDEs (14.321) . we introduce co(t) and Cj{t) through 

Coit) - Coit) = Cjit) - Cj{t) 

= CAit)+Ce{t) + ^cm-{^-o^] 

+ 2N~^ J dxg{x)\^{t,x)\'^ Qn{t,x) , j e N* . 

Thus, co(t) and Cj{t) satisfy 

OO 

9{t) + N^co{t) + ^ n^Cj{t) 

1=1 

(4.33) =-2^ J dxg{x)\^{t,x)\'^gn{t,x) - J dxgix) gn{t, x^ . 

Accordingly, the governing equations of motion are recast into 

idt^{t,x) = [-A +Ve(a:) + g{x) x)\'^ 

(4.34a) +2g{x)gnit, x) - + co(t)] x) , 

idt(j)j{t,x) = [-A + Veix) + 2g{x)^\^{t,x)\‘^ 

+2g{x)gn{t,x) - + Cj{t)] (j)jit,x) 

(4.34b) -A-i$(t,x) J dy ^*{t,y)g{y) . 

Equations (14.331) and (14.341) are the sole consequences of the many-body Schro- 
dinger equation with a state vector of form (14.241) . However, physically it does not 
make sense to have equations of motion for ^j{t^x) with Cj{t) depending on j; each 
^j{t,x) should see the same effective potential. Thus, 9{t) is chosen so that 

each Cj (t) is independent of j . 

We make the simplest choice, Cj{t) = 0 for all j G N. Hence, the global phase 9{t) 
entering (14.241) obeys 

(4.35) 9{t) =-2^ J dxg{x)\^{t,x)\'^gn{t,x) - J dxg{x)gn{t,x)'^ . 

By ()4.34l) . we obtain PDEs (11.11) for ^{t,x) and ^j{t,x) with gs{t,x) = ^|$(t,a;)p, 

(4.36) bj{t) = N-'^J dy^*{t,y)g{y)^\^{t,y)\^^j{t,y) , 

and the thermal cloud density, p„, furnished by (14.101) . For fixed ^ and this 

description concludes our derivation. In section [SJ we propose a closure of the entire 
PDE system via the (equilibrium) Bose-Einstein distribution for 
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5. Stationary states. In this section, we simplify PDEs (EB when the particle 
system occupies stationary states, where corresponds to the lowest macroscopic 
state. Furthermore, we apply the Bose-Einstein distribution for the occupation num¬ 
bers [22] , 

We set 

(5.1) , 

where E is the (lowest) energy per particle of the condensate and Ej is the energy 
per particle of the jth thermally excited state; $ is assumed to be real. By defining 

( 5 . 2 ) ^ = E+^ec, ^iJ = E, + ^ec, 

we re-write PDEs EB as 

(5.3a) {-Aj, -I- Ve{x) + g{x)[gs{x) + 2gn{x)])}^{x) = g.^{x) , 

(5.3b) {-Aa, -I- Ve{x) + 2g{x)[gs{x) + gn{x)]}(j)j{x) - bj $(a;) = gj4>j{x) ; 


(5.4) 


bj = $] = iv ydy 9{y)(t>o{y ); 


gs = and gn = denote the superfiuid and normal fluid densities, and 

C is the static version of (II.2|) . Recall that ||<i>|p = TV and \\(j)j\\ = 1- 

Evidently, fi and fXj are real. Each wavefunction (j)j can be chosen to be real. This 
property follows heuristically by treating $, and hence p„, as fixed in (15.3b)) : Suppose 
(pj is complex. For given functions $, gs and p„, (I5.3bp is satisfied by c{(j)j +4>’j) where 
c is any constant and bj [pj ; $] is replaced by bj = bj [c(((ij -|- ; $]. 

By ansatz (14.241) and formula (I4.35P for the global phase, 0{t), the total energy of 
the Boson system equals 

OO 

£ = e{t) + Y, n^E, (n^'i = N^,E^=E) 

3=0 

OO 

(5.5) = N^E + Y, ^7^3 - 2? 

1=1 

It is plausible to apply the Bose-Einstein distribution, setting 

(5.6) nf = -1)"\ jsN, 

where Uq^ = N^, /3 = l/{kBT) {ksT: Boltzmann energy) and is the 

fugacity of the atomic gas, a Lagrange multiplier (usually denoted by z) [26]. To ease 
notation, we remove the bar and write z instead of z. We impose the constraints [22] 


Axg{x)^{xf‘gn{x) - / Axg{x)gn{xY 


(5.7) TV^ = - 1) \ TV(1-^) - 1) \ 

1=1 

Note that by (15.31) /i and gj can in principle be determined as a function of (C,-^)- 
Hence, EB should yield each of ^ and z as a function of T. 
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6. Periodic microstructure: Two-scale analysis. In this section, we restrict 
attention to a scattering length that varies periodically according to Recall 

that the fast-varying A is 1-periodic. We formally apply two-scale expansion (d to 
PDE system (|5.3p with d by recourse to periodic homogenization [4]. The spatial 
gradient in the governing PDEs is replaced by + (i = xje). We also expand 

OO 

w = '07^°^ -I- e* 777^^^ as e 4- 0 ; zu = ^, ^j, z, bj . 

fc=i 


A two-scale expansion ensues for the densities Qs and Qn- 

In our analysis, we make use of the following standard result pil^. 

Lemma 6.1. Equation —Am = S{-,x), where S{-,x) is 1-periodic, admits a 1- 
periodic solution u(-,x) only if 


(6.1a) 


{S{-,x)) = ( S{x,x) da; = 0 . 


The 1-periodic solution reads as 

(6.1b) M(i, a;) = (—Aj)“^S'(a;, a:)-|-c(x) , 

where c{x) is reasonably arbitrary. 

This lemma, stated here without proof, is a consequence of the Fredholm alter¬ 
native [4]. Any solution of form (I6.1bl) will be referred to as “admissible” |34) . 

6.1. Cascade of equations. Next, we describe the set of equations for coeffi¬ 
cients of two-scale expansion (El- The substitution of expansion (11.41) into PDEs (15.31) 
yields a cascade of equations for and which we write here for fc < 4: 


(6.2a) 

0(e-2) : 

- = 0 =: , 

(6.2b) 

C(e-i) : 

- = 2V, • =: 4'^ : 


0{e°) : 

- A£$(2) = 2V, • - {-A, + Veix) 

(6.2c) 


+ 5o[l + ^(i)](0® + =: , 


0{e) : 

- A£$(3) = 2V, • + Veix) + go[l + A(i)] 

x(p(°)+24°))-^(°)}$W 

(6.2d) 


- { 5 o[l + v4(i)](4) + 2 ^)) - =: 4 ) , 


Oie^) : 

- A£$W = 2 V, • + Veix) + go[l + A(i)] 

x(,(°)+24°))-4)}$(2) 

- { 5 o[l + A(i)](^7« + 24 )) - 

(6.2e) 


- { 5 o[l + A(i)](p( 2 ) + 2 p(f)) - =: 4 ) ; 





BOSE-EINSTEIN CONDENSATION AT FINITE TEMPERATURE 19 


and, for j G N*, 


(6.3a) 

0{e-^) : 

- = 0 > 

(6.3b) 

Cl(e-i) : 

- = 2V, • =: S^'l] , 


0(e°) : 

- + 2V, • ) - {-A, + V,{x) 

(6.3c) 

0(e): 

+ 2go[l -d A(i)](g(°) -d pi°)) - =■ , 

- As</.f = 6f+ 2V, • Vic^f - {-A, + Ve{x) 



+ 2go[i -d A(i)](p(°) -d p1°)) - 

(6.3d) 

0(e2) : 

- {25o[ 1 + ^(x)](eW + 2e(,i)) - /rf =■ ’ 

- Aj<^f ^ = 6f ^$(2) + 6«$(i) + + 2V, • Vi(^f 



- {-A, + V,{x) + 25o[ 1 + A(i)](0i°^ + 1^1°)) - df }<^f ^ 

- {250[1 + ^(x)](5W + el^^) - 

(6.3e) 


- {25o[ 1 + dl(i)](p(^) + g^:'>) - /if =: S^] . 


In the above, gi^\x,x) and git\x,x) are the k-th order coefficients in the expan¬ 
sion for the superfluid and normal fluid densities, viz., 

(6.4a) g’'°\i,x) = , 

(6.4b) x) = x)^-I-a:)$^^^(i, x) , 

(6.4c) ^ ^( 2 )$(o )2 ^( 0 ) ( 2 $( 2 )$(o) $(1)2^ . 


(6.5a) 

(6.5b) 

(6.5c) 


Qn'^ ix,x) = J2 > 

i=i 

OO 

Qn\x,x) = '^[nf'^ 4>f\x,xf + 2nf'’ (t)f\i,x)4>f\x,x)\ , 
i=i 

OO 

ei?> = E{»f <’ + 2d'’ d'’d”’+»f WP'*? +d” 1}. 

i=i 


where are expansion coefficients for Bose-Einstein distribution (15.6^ : 


(6.6a) 

(6.6b) 


(0) 

n) 


= [ 


do)-i /Sm™ _ 


!]■ 


nf ^ = 


( 2 ) 

n) = 




^(0)-lg/3A.f 

[^(o)-v^r _i]2L Vz 


dl) 


do) 


f__ + 3 (2) 

2 .,f0l .^(0) ^ 




(0)-ig/3M™ _ 1 V^: 


( 2 ) 
zW 


dl) 


do) 


(6.6c) 
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Accordingly, constraints (j5.7l) yield the parameters and as a function of 
and For k <2, we obtain the system 


(6.7a) 


(6.7b) 




1 -1 


= ^-E 

i=i 

,,, 




[z(0) - l] 




= -E 




( 0 ) 


- 1 ] 




Ar^(2) = 


^(o) 


'-1 ""W r /^(l)2 (2) 


[z(0)’ 


p/3/^(°) _ 


!]■ 


z(0)^ 


z(o) 


- /3mW 




^(0) 


„( 0 ) 


= -E 


z(0) 


ziO) — 1 V-z 

z 


.( 1 ) 

w 


zW 

z(0) 

2n 




(6.7c) 


^ [z(o)-'e^'^™-l]H 
2\ 

1 J_ 




z(2) a)^ 


^(0)2 z(0) z 


( 1 ) 

M 






2(0)-lg/3M™ _ 1 

Thus, solving for we compute 

(6.8a) z(°) = e^'^™[l + (7V^(°))-i]-' , 

OO 

(6.8b) = Af - ^{[1 + (iV^(o))-i]g^(^r-A‘'°’) _ 1} 

t=i 


■/3m 


( 2 ) 


z(i) 

(6.9a) ^ 

iV^(i) 

(6.9b) 

(o') ciA 

(6.10) +ivcW(z(o)-V'^'">-l)]}+/3^(2)_i(/3^(i))2 . 

Equations (EH) should provide as a function of and {^J■f^}j>l via 

solving a nonlinear equation for Equations ()6.9I) then readily yield (z^), 




= - < 1 + z^'^^e" 






^(0)-lg/3/.f 


1=1 F 


^(0)-lg/3,.f 


(0) 


-1 /3Mf ^i2 




-1]' 


^ [z(o)-ie/5Mf _ 1]' 
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Equation (16.101) can be used in the e-expansion of the second constraint in (15.71) to ob¬ 
tain Finally, are determined from appropriately expanding (EUl), or, 

alternatively, by requiring that ($, =0 be satisfied order by order (see Appendix). 

6.2. Effective equations of motion np to order k = 2. We proceed to 
homogenize PDEs (16.21) and (16.31) in view of (11.31) and (I6.4I) - (I6.10I) under the assump¬ 
tion that A{y) has zero mean, (A) = 0. The resulting equations for the slowly-varying 
parts of and with A: = 0, 1, 2 are displayed in (j6.12|) . (I6.15P and (16.171) . 

By ()6.2al) . (|6.3ap and Lemma 16.11 we trivially have (5'^°^) = 0 = (5'^°]) and 
conclude that 


= /°(x) , = f°{x) , 

by which the zeroth-order densities are gf’\x,x) = Pi{x), for £ = s, n. Similarly, the 
trivially satisfied equations = 0 = yield 

(6.11) $(^^(i,a;) = f{x) , (j)f\x,x) = f}ix) ; 

thus, g^/\x,x) = Pfix). In regard to the density coefficients gf^'^ {k = 0, 1), re¬ 
call (I6.4ap . (|6.4bp . ()6.5ap and (I6.5bp . 

By (I6.2cll and (l6.3cL setting = 0 = (S^j) entails the homogenized equa¬ 

tions of motion for f°{x) and fjix), with ||/°|p = ’n and ||/°|| = 1: 

(6.12a) {-A,, -b Ve{x) + go[g°{x) + 2g°^{x)]}f°{x) = p^°'^f{x) , 

(6.12b) {-A,, -b Ve{x) -b 2go[£>°(x) -b g°^{x)]}f°{x) - bf^f{x) = pf^fj(x) , 

where 6® = goiif^)^, fj) ■ The eigenvalues and of course form part 

of the solution. By Lemma lOI the admissible second-order coefficients of the solution 
read 

(6.13a) <^''^\x,x) = -go[p°{x) + 2pl{x)]f{x)[{-As,)-^A{S:)]+f{x) , 

(6.13b) 4>f\x,x) = -2go[p°(x) -bp°(a;)]/°(x)[(-A£)-iA(i)] -b/j (x) , 

where the (slowly-varying) functions P{x) and ff{x) should be determined. By (I6.4cp 
and (l6.5cL for fc = 2, the second-order coefficients for the densities read 

(6.14a) g^^\x,x) = pI{x) - gopl{x)[{-As,)-^A{i)] , 

(6.14b) gl^'>{x,x) = pI{x) - gofnix)[{-As,)~^A{S;)] , 

where 

pHx) = + 2e^^f\x)fix) + e°H2fix)fix) + f\x)^] , 

p2(x) = 2p°(cr)[p°(x)+2p°(x)] , 

OO 

Plix) = '^{nf^fjixf + 2nfVi (a;)/°(a:) + nf\2ff{x)f°ix) + f-ixf]} , 

1=1 

pUx) = '^[p°{x) -b pI{x)]pI{x) . 
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To find equations of motion for and ff with k = 1, 2, we need to apply 
Lemma [6.11 to the next two higher orders. By (I6.2dl) and (I6.3dl) . the requirement 
{S^^) = 0 = leads to 

{-A^ + Ve{x) + go[p°s{x) + 2p°(a;)] - 
(6.15a) = - go[plix) + 2pl^{x)]}f°{x) , 

{-A^ + Ve{x) + 2gQ[p°{x) + pI{x)] - pf^}f}{x) - f 

(6.15b) = - 2g^[pl{x) + pi{x)\}f°{x) , 

where and are subject to (/^,/°) = 0 and {fj,fj) = 0; and b^^^ is found 
from imposing + (/°, fj) = 0 (see Appendix for a derivation): 

= {(/\p°/°) + {fj,P°f°) + if^plfj)} 

= iv-l5o{e^°)((/°)^//) + 3e(°)(/^(/°)V°)+e^'^((/°)^/°)} • 

By Lemma l6. 11 the third-order coefficients in the e-expansion of d> and (j)j read 

<i>(3)(i,a.) = - 25 oV 4 (p°(x) +2p0(x))/°(x)] • y,{AfA) 

- 9o[{p°ix)+‘2Pnix))fix) + {plix) +2pl{x))f{x)] 

(6.16a) X {(-As)”^^} -H f{x) , 

(l)f\x,x) = -45oV^[(p°(x) -hp°(x))/°(x)] - 

- 25 o[(p°(a:) -Hp°(x))//(x) -F {pl(x) + pl{x)) f°ix)] 

(6.16b) X {(-Aj)"^A} -P ff{x) . 

We continue our homogenization program, proceeding to the next higher order. 
By (I6.2ep and (|6.3ep . setting (5'^'^^) = 0 = (*5'^'^]) yields the following PDEs. 

{-A^ -P Ve(x) + go[p°s{x) + 2p°(x)] - p^^^fix) 

= { - 9o{pI{x) + 2pl{x)) +p^^^}f{x) 

+ {9l\\M-i{Psix) + 2p°(x))^ - 9l\\M-i{plix) + 2pl{x)) 

(6.17a) - go{plix)+2pl{x)) + p^‘^^}f{x) , 

{-A, + Veix) + 2go[p°(x) + (x)] - pfjffix) - bf f - bf^f^ - bf f 

= { - ‘^ 9 q{p\{x) + p\{x)) A pf'>)f]{x) 

+ {'^9l\\A\W{p'l(.x) + Pn{x)f - 2gl\\A\W{pl{x) -H p^(x)) 

(6.17b) - 2go{pl{x) A pl{x)) -h /rf^}/°(x) , 

where ||A||_i is the iJ“^-norm of the 1-periodic, fast-varying function A(x); recall 
that (A) = 0. The factor b^p is determined directly from (15.41) which is consistent 
with imposing the condition (/°, /J) -|- + {P,fj) = 0; see (IA.5I) in Appendix. 

The parameters and p^ are subject to (/°,d’(^)) -|- (/^,/^) -I- /°) = 0 
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and = 0 which, by Lemma lA.ll of Appendix, yield 

||/i||2 2(/°, p) = 0 and ||//p + 2(/°, /|) = 0, respectively. 

The main results of our formal procedure are summarized as follows. 

Proposition 6.2. Consider governing equations (lOl) for $(a::) and 4>j{x) (j € 
N* j under spatially-varying scattering length (131). Two-scale expansion dUD entails 
x) = P{x) and = /^(x), which depend only on the slow variable, 

while and (f>^^\x,x) are given by (16.131) . The slow-varying functions f^{x) 

and fj{x) (k = 0, 1, 2) satisfy (16.121) . (16.151) and (16.171) . In these equations, the pa¬ 
rameters and appearing in the density coefficients and pit'' are functions 
and pty' according to (16.81) - (16.101) with (15.71) . 

6.3. Energy of Boson system. We conclude the homogenization program by 
discussing the total energy (1^ . The requisite energy parameters are E and Ej-, 
cf. section [5l By (15.2|) . we deduce that 

(6.18) A = + ... , Ej = E^°^ + eE^^'> + e'^Ef^ + ... , 

where 

E(0) = po) _ l^(0)2^(0) ^ 

eW = p(i) - i [^(0)^(1) + , 

(6.19) e( 2) = p(2) _ i [^(0)^(2) + ; 

and similarly for Ej via replacement of by In the above, result from 

the e-expansion of the static version of m via Lemma |A.2l of Appendix; viz.. 


^ = ^(0) + g^(l) + g2^(2) ^ ^ ^ _ 

= N-^9o{\mY + Mf\ iff) + [4(/^ im + 6||/V°f 
-4go\\A\\\\\PpTMinY]} • 


Consequently, by (OD the total energy reads £ = -\-eP'-'^ -\-e ^-\-..., where 

£:(0) ^ ^^(0) ^ ^ ^(0)^(0) _ ^0^ _ ^^11^0 ||2 ^ 

3 

£:(1) = + Ef^ + nf''Ef^) 

3 

-‘^9o{{pIijP°s) + (pl^Pn) + (PnfPli)) ) 

£:(2) = E^^^ -b -b 

3 


2go{(p°,p^) + {pI^pD + {p°,pI) + {pl,Pn) 
■ 9o\\Ap{{Pn,fs) + (Ps.Pn) + (P°,Pn))} ■ 




( 6 . 20 ) 
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7. Discussion and conclusion. In this article, we introduced a formalism for 
Bose-Einstein condensates in an inhomogeneous scattering environment at sufficiently 
low, finite temperatures. By starting with the Hamiltonian description of N repul¬ 
sively interacting Bosons forming a dilute gas in a trap, we formally derived mean 
field evolution equations for the condensate wave function and the single-particle wave 
functions of thermally excited states. Our main assumption is that the temperature 
is finite but lies sufficiently below the phase transition point. In our model, the 
scattering length, which expresses the strength of the pairwise particle interactions, 
is phenomenologically included in the Hamiltonian as a spatially varying, periodic 
function of subscale e. 

First, we reduced the many-body problem to lower-dimensional evolution PDEs 
on the basis of a perturbative treatment of the many-body Hamiltonian in Fock space. 
This methodology is distinctly different from variational approaches for stationary 
settings, e.g. [IHHH]; and enabled us to derive mean field laws for time-dependent 
one-particle wave functions. Our analysis is a direct extension to finite temperatures 
of techniques applied to the case with zero temperature |23ll46j . In the present work, 
we restricted attention to formal derivations, by invoking uncontrolled ansatz (14.241) 
for the A^-body Schrodinger state vector. The particle occupation numbers, nj, for 
thermally excited states are treated as parameters in the evolution laws. 

Second, we homogenized the resulting PDEs up to second order in e for stationary 
states. This procedure relies on a scale separation for the condensate wave function 
and the one-particle wave functions of the thermally excited states. As a result, we 
obtained a system of e-independent equations of motion for the slowly-varying parts 
of the single-particle wave functions. In these effective equations, the oscillatory zero- 
mean ingredient. A, of the scattering length contributes terms proportional to its 
norm squared, analogously to the case with zero temperature [34) . 

There is a number of issues that remain unresolved within our treatment. In opti¬ 
cal traps, the rapid spatial variation of the scattering length is driven by the external 
potential, Ve(x), set by appropriate laser fields [9]. This suggests that, in realistic 
settings, this Ve(x) should also vary spatially with the same subscale, e. The joint 
effect of a rapidly varying 14 and the periodic scattering length, a, is the subject of 
work in progress. The current treatment, placing emphasis on the extraction of low¬ 
dimensional laws of motion, has not addressed the issue of approximate solutions to 
these equations, e.g., in the Thomas-Fermi regime [10]; in this regime, it is possible to 
derive an approximate description for the superfluid density combined with an approx¬ 
imation for the normal fluid density. This task is left for future work. The spatially 
dependent scattering length is introduced as an ad hoc parameter in the microscopic 
Hamiltonian. The emergence of such a parameter from the limiting procedure lies 
beyond our present scope. The finite-temperature approach here introduces the oc¬ 
cupation numbers, nj, of one-particle thermally excited states as given parameters in 
the quantum dynamics; the derivation of these parameters was not addressed. The 
present formalism leaves out two effects: quantum depletion of the condensate due to 
pair excitation, which calls for modifying ansatz (14.241) : and the effect of dissipation 
due to couplings of the condensate with the environment. The effect of quantum de¬ 
pletion, brought about by the operator V 21 of (14.131) . should be tractable through the 
introduction of the suitable pair-excitation kernel into the ansatz for the many-body 
Schrodinger state vector [46]. This approach, which transcends the usual mean held 
limit, is a promising direction of research. Finally, we believe that the approxima¬ 
tions for the A^-body Hamiltonian are valid over some moderate time scale; the precise 
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characterization of this time regime has been left unresolved. 

(k) 

Appendix A. On the computation of o'^ , for fc = 1, 2. In this appendix, 

we compute the coefficients in the e-expansion of formula (15.4p with fc = 1, 2. For 
this purpose, we need a few results pertaining to classical asymptotics of integrals of 
highly oscillatory functions [33]. 

Lemma A.l. Consider the integrahle function ip : —>■ M and the bounded, 

1-periodic P : > R with (P) = 0 (d > 1). Suppose that ip has m summable 

derivatives, to £ N*, which vanish at infinity; then, 

(A.l) [ p( —] p{x) dx = 0{e^) as e I 0 . 

Js.'‘ V e / 

In the special case with smooth p (p G 

(A. 2 ) 


for any k G N. 

The proof of this lemma can be established via successive integrations by parts [34] , 
and is omitted here. 

More generally, for 1-periodic and bounded P : ^ K with (P) 0, one 

can write P{y) = (P) -I- P{y), where P{y) satisfies the hypotheses of Lemma [A.11 
particularly, (P) = 0. Thus, we reach the following conclusion. 

Lemma A. 2. Consider the integrable p : ^ R. and the bounded, 1-periodic 

P : —>■ R. with finite (P), (P) ^0 (d > 1). Suppose that p has TO summable 

derivatives, to £ N*, which vanish at infinity; then, 

(A.3) f P(—j p{x) dx = (P) f p{x) dx-\-0{e^) as e I 0 . 


We proceed to describe the computation of (fc = 1, 2). By inspection of (15.4p 

along with the formulas for and provided in section 16.21 we need to use 
Lemma lA.ll for the bounded (and smooth) P{y) = A{y) and Lemma IA.2I for the 
bounded P{y) = A( 2 /)[(—Ay)“^A(?/)], y = x/e. On the other hand, p[x) contains 
products of /^’s and /j^’s, each of which is assumed to be smooth. Note that 

(A[(-A)-M]) = ((-A)-iA, A) = ||Af_i , 

which is the P“^-norm squared of the (1-periodic, zero-mean) A(y). 

Without further ado, by $ ^ /° -|- ef^ and ^ fj efj (15.4p 

yields 

= A-igo lim / [1 + A{x/e)] [f}{x)f{xf + 3 f {x)^ f\x) f° (x)] 

+ CV°(a^)V°(a^)} 

= [((/°)^ fj) + 3((/°)V°, /')] + e^'^((/°)^ fj)} 


(A.4) 
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and 



|(^(0)+,^(l)+,2^(2)) f [l + A{x/e)] 

I JR3 


X [f{x) + ef\x) + e^^^^\x/e,x)f 


X 



(A.5) 


[((/°)^ ff) + 3(((/°)^ /Vj) + ((/')^ /V°) 

+ ((/°)V°, r)) - 5oPf (/°, (5p° + 8p°)(/°)3)] 
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